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ABSTRACT 



Context. Mapping the brightness distribution of exoplanets is the next frontier for exoplanet infrared photometry studies. For tidally- 
locked hot Jupiters that transit and are eclipsed by their host star with non-zero impact parameter, the first steps are now possible. 
Aims. The aim is to use eclipse scanning from occultation ingress/egress and phase curve measurements to constrain exoplanet large- 
scale brightness structure. 

Methods. We use archived Spitzer/IRAC 8yL/m data of HD 189733 in a global MCMC procedure encompassing six transits, eight 
secondary eclipses, and a phase curve in a two-step analysis. The first step derives the planet-star system parameters. The second step 
investigates the structure found in eclipse scanning, using the previous planet-star system parameter derivation as Gaussian priors. 
Results. We find a 5-cr deviation from the expected occultation ingress/egress shape for a uniform brightness disk, and demonstrate 
that this is dominated by large-scale brightness structure and not an occultation timing off'set due to a non-zero eccentricity. Our 
analysis yields a 2D brightness temperature distribution showing a large-scale asymmetric hot spot whose finer structure is limited by 
the data quality and planet orbit geometry. We also present an improved upper limit for eccentricity, e < 0.0081 (95% confidence). 
Conclusions. Reanalysis of archived HD 189733 data revealed brightness structure by using global analysis that mitigated systemat- 
ics. Future eclipse scanning observations of the same exoplanet at other wavelengths will probe diff'erent atmosphere layers, ultimately 
generating a large-scale 3D map. 

Key words, eclipses - planets and satellites : individual : HD 189733b - techniques : (photometric) 



1. Introduction 

Among the hundreds exoplanets detected to dat^ it is the subset 
of transiting exoplanets that can be most extensively character- 
ized with current technology. Beyond derivation of the exoplanet 
orbital inclination and density (e.g., Winn 2010 ), transiting plan- 
ets are key objects because their atmospheres are observationally 
accessible through transit transmission and occultation emission 
spectrophotometry (see e.g., Seager & Deming 2010, and refer- 
ences therein). In addition, exoplanet phase-dependent infrared 
flux modulations provide observational constraints on the ex- 
oplanet atmosphere p roperties; from phase curve IR measure- 
ments, Knut son et al.|p007 ) derived the first longitudinal bright- 
ness distribution of an exoplanet. This longitudinal map indi- 
cated a hot spot within the exoplanet atmosphere, in particu- 
lar an offset hot spot in agreement with predicted super-rotating 
equatorial jet for hot Jupiters (e.g., [Showman & Guillot||2002| ). 
Related, significant theoretical developments have also been 
achieved in modeling exoplanet atmospheric hydrodynamic flow 
with radiative transfer (e.g., [Showman et al.|[2Q09t|Rauscher &| 
|Menou|2010||Dobbs-Dixon et al.|201Qt [Moses et al.|2Qll| ). 

The observation of localized, specific spatial features, such 
as hot spots or cold vortices, is essential for constraining atmo- 



^ For an up-to-date list see the Extrasolar Planets Encyclopaedia: 
http : //exoplanet . eu| 



spheric structure that arises from complex hydrodynamic pro- 
cesses. Eclipses have proved to be powerful tools for "spatially 
resolving" distant objects, including binary stars (e.g., Warner] 
et al. 197 1) and accretion disks (e.g.. Home 1985) . Previous the- 
oretical studies introduced the potential of eclipse scannin^see 
Figure[T]) for exoplanets in order to disentangle atmos pheric cir- 
culation regimes (e.g., Williams et al. 2006; Rauscher et al 



2007 ). In the first direct mention of eclipse scanning, Agol et al 



(,2Q10j dete cted the "uniform time oflTset" (defined by Williams 



etal. 2006 see our Subsection 2.2 ) expected for the HD 189733b 



hot spot based on the longitudinal map derived by |Knutson et al. 
12007)). 

Ideally a tidally-locked transiting and occulted exoplanet 
with a non-zero impact parameter can enable a coarse 2D surface 
brightness map. As represented in Figure[T] such an exoplanet is 
scanned through several processes along its orbit. During occul- 
tation ingress/egress, the planet is gradually masked/unmasked 
by its host star. In addition, one may derive the flux in longitu- 
dinal slices from the phase curve as the planet revolves around 
the star, as long as the exoplanet is tidally locked. The three dif- 
ferent scanning components (ingress, egress, phase curve) pro- 
vide us with complementary pieces of information that could ul- 



^ Eclipse scanning is the process by which a body gradually masks 
another body. 
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Fig. 1. Schematic description of the different scanning processes observable for a occulted exoplanet. 



timately lead to determine the brightness distribution over a spe- 
cific "grid" (e.g., see the component labeled ''totaF in Figure[T]). 

Exoplanet thermal infrared observations aiming at character- 
izing planetary brightness distribution are growing. Currently, 
Spitzer has observed thermal phases curves for eight diff'erent 
exoplanets as well as the infrared occultations of over two dozen 
exoplanets. Among these, HD 189733b ( |Bouchy et aLl|2Q05| ) is 
arguably the most favorable transiting planet for detailed obser- 
vational atmosphere studies: its K star is the closest to Earth 
with a transiting hot Jupiter. This means the star is bright and 
the eclipses are relatively deep yielding favorable SNR datasets. 
As such, HD 189733b represents a "Rosetta Stone" for the field 
of exopl anetology with: one of the highest SNR sec ondary 
eclipses ( Deming et al.||2006[ Charbonneau et^|2QQ8 ); phase 
curve observations ( [Knutson et al.||20Q7[ |2Q09| ); numerous at- 



mospheric observational characterizations (e. g.,|Grillmair et aL 
2007; Tinetti et al. 2007; Redfield et al. 2008llS wain etal.|20081_ 
Madhusudhan & Seager 2009; Desert et al. 2009; ?; |Deroo etal.| 
|2010). Although HD 189733b atmospheric models are in quali- 
tative agreement with observations, important discrepancies re- 
main between simulated and observed light curves as well as 
emission spectra (see e.g.. Showman et al.|2009| ). 

We undertake a consistent and global analysis of all 
HD 189733b public photometry obtained with the Spitzer Space 
Telescope ( Werner et al. 2004 ) to assess the validity of published 
inferences and detect effects that would remain unnoticed in the 
frame of individual data sets analysis (de Wit et al., in prep). 
In this paper, we present the first eclipse scanning with a detec- 
tion of an occultation ingress/egress deviation from a uniform 
brightness disk at the S-cr level. However, the brightness deter- 
mination over a specific grid (see Figure[T]) is not yet obtainable 
because 1) the data is not of high quality enough: there will be 
no significance on a grid cell and 2) to make a 2D map a higher 
level of precision of the system parameters is needed. In par- 
ticular for the eccentricity that is highly degenerated with the 
brightness distribution. For that reason we proposed another pro- 
cedure to constrain consistently the eccentricity-brightness de- 
generacy. The main point of the paper is our new procedure on 
how to handle the eccentricity-brightness distribution degener- 
acy and what features in a 2D surface brightness distribution can 
be identified. In particular, we simultaneously derive the domi- 
nant two-dimensional pattern of HD 189733b's atmosphere and 
a highly precise upper limit on its orbital eccentricity, by com- 
bining the occulta tion ingress/egress dev iation detection with the 
phase-curve from [Knutson et al.](|2007|). 



At the time of submission, we learned about a similar work 
by |Majeau et al. ( 2012j ), hereafter Ml 2, focusing on the deriva- 
tion of HD 189733b 2D eclipse map, using the same data set 
but a different data analysis. In contrast to Ml 2 who present a 
2D map, we focus on the dominant 2D surface brightness pat- 
terns, because finding a unique 2D map requires higher qual- 
ity data than available. Furthermore, we present a 5-cr deviation 
from the expected occultation ingress/egress shape for a uni- 
form brightness disk, a critical component in deriving the 2D 
dominant brightness pattern. Our work differs from Ml 2 in two 
ways. First, a significant point is that the planet surface bright- 
ness distribution is degenerate with planet orbital eccentricity 
which is not measured. Our work highly constrains this degen- 
eracy. Second, and related to the first point, we propagate the 
system parameter uncertainty as it affects the precision to which 
the 2D surface brightness patterns can be derived. We further 
discuss a comparison of the techniques and results in Sec.|4] 

We begin with a su mma ry description of the 8yum Spitzer 
data sets. In Subsection 2J_ we present our photometric data 
analysis and the occultation ingress/egress deviation from a uni- 
form brightness disk. In Subsection 2.2 we investigate the possi- 
ble origins of the ingress/egress structure and describe our proce- 
dure for consistently disentangleing the causes. We present our 
results in Section|3j including the HD 189733b brightness dis- 
tribution at 8yum and a precise eccentricity upper limit. Then 
we discuss the robustness of our detection and the necessity of 
global analysis to mitigate the impacts of systematics, to both en- 
able detection of previously unnoticed effects and constrain de- 
generacies. We conclude by discussing the implications for mea- 
suring the dominant 3D brightness features: each wavelength 
corresponds to a different atmospheric layer in altitude and so 
combining 2D results at different wavelengths can lead to an ex- 
traction of 3D map features, highlighting tremendous potential 
for resolving distant worlds. 



2. Data Analysis 

We divide our work in two major steps. First, we perform a 
global determination of the system parameters (the orbital and 
physical parameters of the star and planet) based on the IRAC 8 
yum transit and occultation photometry. The parameter determi- 
nation lead to our detection of anomalous ingress/egress shapes 
that are presented at the end of Subsection |2.1| Using the sys- 
tem parameter estimates as Gaussian priors, we simultaneously 
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analyze the Syum Spitzer IRAC photometry of HD 189733, in- 
cluding occultation and phase curve photometry. 

2.1. System Parameters Estimation 
2.1 .1 . Data description and reduction 

The eight secondary eclipses and the six transits of HD 189733b 
used in our system parameter estimation are listed in Table[2] 
The data was obtained from November 2005 to June 2008 with 



the Infrared Array Camera (IRAC: [Fazio et al.]|2004| ) at 8yum. 
The dat^ used in this study, calibrated by the Spitzer pipeline 
version SI 8. 18.0, are described (by Astronomical Observation 
Requests, hereafter AOR) in Table|2] The new SI 8. 18.0 version 
enables improvements in the quality of data reduction over the 
original published data sets which used older Spitzer pipeline 
versionqj 

The data reduction follows conventional procedures (e.g., 
[Gillon et al.||2010b[ hereafter GIO, and references therein) and 
is performed individually and consistently for each AOR (de 
Wit et al., in prep). For each AOR, we first convert fluxes from 
the Spitzer units of specific intensity (MJy/sr) to photon counts, 
then perform aperture photometry on each image with the IRAl|j 
/DAGPHOT software (Stetson|[T9S7j). We estimate the best aper- 
ture radius (see Table[2]) based on the instrument point-spread 
function (PSiQ), HD189733b's, its host star and the background 
contributions. We estimate the PSF center by fitting a Gaussian 
profile to each image. A mean sky background is measured in 
an annulus extending from 10 to 16 pixels from the PSF cen- 
ter, and subtracted from the measured flux of each frame. We 
discard the first 30 minutes of the data to allow detector stabi- 
lization. According to the x-y distribution of the PSF centers we 
reject the few significant outliers to the bulk of the distribution. 
For each set of 64 subarray images, we discard measurements 
with discrepant values of flux, background and PSF center posi- 
tions using a lOcr median clipping. The resulting flux values are 
averaged across each set of 64 subarray images and we take the 
photometric error as the error on the average flux measurement. 

2.1.2. Photometry data analysis 

We now describe the procedure we use to analyze the eight oc- 
cultations and six transits obtained in IRAC 8 jim channel in or- 
der to constrain HD 189733b system parameters (see Table[2]). 
We emphasize that we simultaneously analyze the whole data 
set in what we term a "global analysis", instead of combining 
the each separately analyzed transit or eclipse events. The global 
analysis approach helps to mitigate the impacts of systematics 
and therefore enables detection of previously unnoticed efl'ects 
in the data. For this first part of the analysis, we model the exo- 
planet as a disk of uniform brightness distribution. 

We use an adaptive Markov Chain Monte Carlo (MCMC; see 
e.g. |Gregory|2005t [Ford|2006| ) algorithm. MCMC is a Bayesian 
inference method based on stochastic simulations that sample 



^ Data available in the form of Basic Calibrated Data (BCD) on 
the Infrared Science Archive : |http : //sha . ipac . caltech . edu/| 
iap plicatio ns/Spitzer/SHA// | 



http:// irsa . ipac . caltech. edu/data/SPITZER/docs/ 



'irac/iracinstrumenthandbook/73 

^ IRAF is distributed by the National Optical Astronomy 
Observatory, which is operated by the Association of Universities 
for Research in Astronomy, Inc., under cooperative agreement with the 
National Science Foundation.^ 

^ |http : //ssc . spitzer . caltech . edu/irac/psf . html 



the posterior probability distribution of adjusted parameters for 
a given model. We use here the implementation presented in de- 
tail in Gillon et al.] ( [2009 , 2010a b). More specifically, this im- 
plementation uses Kepleri an orbits and models th e eclipse pho- 
tometry using the model ofl Mandel & Agol| ( |2002| ) multiplied by 
a baseline, diff'erent for each time-series, to take into account the 
systematics (see below). 



Eclipse model & limb-darkening 

We model the secondary eclipse assuming the planet to be a uni- 
form disk, and the transit assuming a quadratic limb-darkening 
law. We draw the limb-darkening coefl&cients from the theoret- 
ical tables of Claret & Bloem en| ([20TT]), ui = 0.0467 ± 0.0009 
and U2 = 0.0979 ±0.0015, based on the spectroscopic parame- 
= 5050 ± 50^ log g = 4.61± 0.03 and 
and references therein). 



ters of HD 189733 (Tefi 

[ f ] - 0.03 ± 0.05, see|Southworth 



2010 



We add this a priori knowledge, using as additional jump param- 
eter^ the combinations ci = 2ui + U2 and C2 = ui - 2u2, as a 
Bayesian penalty to our merit function as described in GIO. The 
impact of this coefficient choice does not aff'ect our results (see 
Sec.[3]). 



Systematics 

IRAC instrumental systematics variation of the observed flux, 
such as pixel-phas e or detector ramp, are well-documented (e.g.. 
Desert et al.|2009[ and references therein). At 8 yum. Si: As-based 
detector showed a uniform intrapixel sensitivity (i.e. negligible 
pixel-phase eff'ect) but te mporal evolution of their p ixels gain 
(detector ramp). Following [Charbonneau et al^p008| ) we model 
locally the detector ramp as a quadratic function of ln(dt). In ad- 
dition, we take into account known low-frequency noise sources 
(instrumental and stellar) with a second order time-dependent 
polynomial. 

We use linear baseline models and determine their coef- 
ficients by linear least squares minimization at each step of 
the MCMC. For this purpose, we applied the S ingular Value 
Decomposition (SVD) method ( [Press et al.|1992 ). 



Correlated noise 

To obtain reliable error bars on our parameters, we take i nto ac- 
count the correlated noise following a procedure similar to | Winn 



et al. (2008]). For each light curve, we estimated the standard de- 



viation of the best-fitting solution residuals for time bins ranging 
from 3.5 to 30 minutes. For each time bin, the following factor 
jSred is determined 



fired - 



(Tn N(M - 1) 



M 



(1) 



where N is the mean number of points in each bin, M is 
the number of bins, and cri and cr^ are respectively the standard 
deviation of the unbinned and binned residuals. The largest value 
obtained with the different time bins is used to multiply the error 
bars of the measurements. 



^ Jump parameters are the model parameters that are randomly per- 
turbed at each step of the MCMC method 
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Table 1. Exoplanet and star parameters 



Parameters (units) 


Median + 1 cr 




0.023985±g:S 






P(days) 


2.2185740±g:— 


W(days) 


07513+^-^^^^^ 


Toihjd - 24539S0) 


8 SH'^OSI _i_0.000066 


FplFi.\%Hm 




a(AU) 


03143+^-^^^^^ 


in 


79 ,0.058 


P*(Po) 


1 905+0037 

^•^^-^-0.024 




762+^-^^^ 


Rp(Rj) 


1 152+^-^2^ 

^•^^^-0.028 



Jump parameters & priors 

In our analysis, we used the following jump parameters; the 
planet/star area ratio (Rp/Ri,)^, the orbital period P, the transit 
duration (from the first to last contact) W, the impact parameter 
b = a/Ri, cos /, cos oj and ^fe sin oj and the time of minimum 
light Tq. We assumed a uniform prior distribution for all these 
jump parameters and draw at each step a random based on 
the Gaussian prior: = 0.84 ± Q.Q6M0 ( [South worth|201Q| ). In 
this first analysis, we also constrained the eccentricity to zero 
(see below and Sec. [3]). 

2.1.3. Intermediate results 

We present the system parameter estimation results here, as they 
are the starting point for the second analysis step focussing on 
HD 189733b eclipse scanning (see Subsection|2.2[). 



System parameters 

We present the system parameters resulting from our global 
MCMC simulation in Table[T] by the median values and 68.3% 
probability interval for the jump and physical parameters. We 
estimate these parameters similarly to |Agol et a l.|(|201Q]), setting 
e = based on the small inferred value of e cos co and theo- 
retical predictions (see Sec. [3]). We further test the incidence of 
the system parameter estimation based on the ^ = assumption: 
the conclusions presented in Section|3]are similar when based on 
system parameters estimated with e freely varying or with e set 
to zero. We compute our best-fit parameters based on a global 
MCMC simulation of the ten AORs (see Table|2]), in opposition 
to Agol et al. ( 2Q10| ) that took these as a weighted mean of in- 



dividual transit or eclipse analysis. Although our global analysis 
is more consistent to extract the information of each AOR and 
thus robust to systematics, our best-fit parameters are in good 
agreement with | Agol et al.| ( [20TQl ). 



Eclipse scanning 

We have detected structure in the HD 189733b occultation 
ingress/egress light curves that deviates at the 5cr level from the 
expected Ingres s/egress shape of a gradually occulted uniformly 
bright disk. We determine the significance of this structure as 



jSred (here on average of 1.3) to take into account the correlated 
noise eff'ects on our detection. Figure[2] shows the IRAC pho- 
tometry of the eclipses ingress/egress corrected for the system- 
atics and binned per 1 minute, with the best-fitting eclipse model 
superimposed. Our detection occultation ingress/egress devia- 
tion detection is in agr eement with the expected signature of the 
HD 189733b hot spot ([Knutson et 'aLl[2Q07 ). We explore other 



possible origins for this structure in Subsection 2.2 



2.2. Eclipse scanning analysis 

The second part of the analysis aims on investigating the ori- 
gin of the occultation ingress/egress deviation from a uniform 
brightness disk, i.e. the 5-cr detection of structure in the occulta- 
tion ingress/egress. 



2.2.1 . Eclipse scanning: possible origins 

The primary potential origins of the occultation ingress/egress 
structure detected are: (1) non- spherical exoplanet; (2) small, 
but non-zero, eccentricity; (3) non-uniform brightness distri- 
bution. (1) We reject the non-spherical exoplanet concept be- 
cause we have detected no significant structure in transit 
ingress/egress (Figure[2j in particular bottom-right panel). The 
spherical planet finding is in ag reement with cu rrent constraints 
on the HD 189733b oblateness ( Carter & Winn|20 1 0) and wind- 
driven shape (expected to introduce light-curve deviation below 
lOppm, see |Bames et al.|2009| ). For HD 189733b shape-induced 
structure are expected to be about one order of magnitude larger 
in primary than in secondary eclipse (Ip/h, Ip,h respectively 
the exoplanet and the star mean intensity at 8 yum). We therefore 
take the planet to be spherical. (2 and 3) A non-zero eccentricity 
can introduce an apparent "timing off'set" to the occultation of 
a uniform disk. This can mimic a non-uniform brightness dis- 
tribution that is also expressed by structure in the occultation 
ingress/egress (see William s et al.|2006| ). In particular for nearly 
circular orbits, o ccultations o ccur Ate ~ f [1 + ^^coso;] after 
transits (see e.g., |Winn|2010 ), leading to a biased estimation of 
^/e cos oj and ^/e sin oj due to apparent timing off'set (see Sec. [3]). 

Before turning to our occultation mapping we recall that ex- 
oplanet scanning involves three diff'erent components: ingress, 
egress and phase curve. They provide complementary pieces 
of information (for a tidally-locked planet with non-zero im- 
pact parameter) that could ultimately lead to a determination of 
the brightness distribution over a specific "grid" (see Figure[T]). 
However a 2D map is not yet obtainable because 1) the data is 
not of high enough quality: there will be no significance in sig- 
nal on a grid cell and 2) to make a 2D map a higher level of 
precision of system parameters is needed, especially for the ec- 
centricity, than is available. Here we provide a careful analysis of 
the brightness distribution-eccentricity degeneracy and use that 
to extract dominant features in a 2D brightness distribution. The 
analysis also leads to confirmation that the ingress/egress bright- 
ness anomaly is real and not solely due to a timing off'set from 
eccentricity. 



I^ieingress ^ " I^ieegress whcrc 7, and o", are the measure- 
ment and its error bar at time /. The error bars are rescaled by 



2.2.2. Data analysis using a non-uniform brightness source 
occultation model 

We turn to describe the procedure we used to constrain the 
eccentricity -brightness degeneracy using the complementary in- 
formation from HD 189733b occultation and phase curve mea- 
surements. The data we use in this second step of the analysis are 



4 




Fig. 2. IRAC 8 fim HD 189733b occultation and transit ingress/egress. Left: Folded occultations Ingres s/egress deviates from the 
uniform brightness disk eclipse model, highlighting the eclipse scanning of HD 189733b dayside. Right: Folded transits showing 
no significant structure during ingress/egress (i.e. no significant deviation to the transit of a perfect sphere). Top: Global eclipse 
photometry binned per 1 minute and corrected for the systematics with the best-fitting eclipse model superimposed (in red). Bottom: 
Combined residuals. 



the detrended and folded occultations obtained from the first step 
of our analysis (from S ub section ] 2. 1| ), and the phase curve from 
( Knutson et al.||2QQ7| ) corrected for stellar variability (see |Agol 
et al.|201Q ). We discard the first third of the phase curve since it 
is strongly affected by systematics, incl. the ramp correction. 

The procedure we use to constrain the eccentricity- 
brightness degeneracy uses a second MCMC implementation 
that diff'ers from the one introduced in Subsection l2.1l in that it 
uses non-uniform brightness distribution for the exoplanet. An 
additional major diff'erence is that we use as additional input the 
planet- star system parameter estimation from the first step of our 
analysis (see Subsection 2J_ -h Table[T]). The second MCMC al- 
gorithm models the light curve photometry by taking a bright- 
ness distribution based on simple analytic brightness models (see 
below) and then integrates the exoplanet flux. Approximations 
in the model include: ignoring the time variability of the target 
atmosphere (in line with atmospheric models, e.g., [Cooper & 
Showmaii||20Q5 ); ignoring the planet limb darkening (expected 
to be weak for hot Jupiters, in particular at 8 yum); and assum- 
ing HD 189733b to be tidally l ocked (as non-e ccentric old hot 
Jupiters are expected to be, e.g., Fabry cky|2Q 10 ). 

We perform the photometry by sampling the exoplanet sur- 
face with a grid of 2N points in longitude (0) and N points in 
latitude (6). We fixN = 100 to mitigate numerical eff'ects. To es- 
timate the broad patterns of HD 189733b brightness distribution, 
we use toy models as; 



Ti((P,0) = /iCexp[-0?]cos^^o+/o, 

h cos ^(po cos ^6o + /o if(po > 
Ii cos^0o cos^^o + Iq if(po < 



r2(0,^) 



(2) 
(3) 



r m - exp[-(^)2] exp[-(^)2] + k ifcp. > 
- I /, exp[-(|)^] exp[-(^)^] + /o if CP. < ' 

where cpo and 6o are respectively the longitude and lati- 
tude relative to the position of the model maximum, i.e. cpo = 
f((f),A(p,Model,a,p,y) and Oo = f(e, AO, Model,a,j3,y). Acp 
and are respectively the longitudinal and the latitudinal shift 



of the peak model from the substellar point. a^jS and y constrain 
the model extension. Our toy models add to traditional MCMC 
algorithm five jump parameters (A(p, A6,a,j3,y) and two linear 
coefficients (/i, /o), respectively to the amplitudes of a "hot/cold 
spot" mode (r/(0, 6) - Iq) and the constant mode (/q). 

At each step of the MCMC, we first model the target orbits 
and the model light curves. We then derive the mode amplitudes 
using the SVD method. 

Note that we use three separate models in order to mitigate 
the model-dependence of our results. As discussed in Sec. [3] 
each of the models lead to similar conclusions regarding the con- 
straint on the eccentricity-brightness degeneracy and the major 
patterns of the HD 189733b brightness distribution. 



3. Results 

Non-uniform brightness from occultation ingress/egress devi- 
ation. 

Our first result is a robust, 5-cr detection of the occultation 
ingress/egress deviation from the eclipse scanning of a uni- 
formly bright disk. This result was introduced in Subsection |2.1[ 
Our upper limit on eccentricity combined with the spherical 
planet limit (based on a uniform transit ingress/egress), means 
that the occultation ingress/egress anomaly is mainly due to a 
non-uniform brightness distribution. On it's own the anomaly 
detected in occultation ingress and egress separately provide an 
average brightness distribution in slices on the planet as shown 
in Figure[T] however it does not constrain the brightness distri- 
bution in a 2D fashion. 

We test the robustness of the occultation ingress/egress 
anomaly detection against various eff'ects including the baseline 
models and the limb darkening. In particular we also validate the 
result's independence to subsets of AORs by analyzing diff'erent 
subsets of the eight eclipses. The diff'erent MCMC simulations 
show that there is no eff'ect to within 1 cr. 
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Dominant features in a 2D brightness distribution at 8yum 

We have extracted a major spatial feature in HD 189733b bright- 
ness distribution, that corresponds to a hot spot. There is a double 
asymmetry in North-South (17 ± 10°) and East- West (24 ± 11°) 
within the brightness distribution. Figure[3] left panel, shows the 
extracted HD 189733b large-scale relative brightness (i.e. over 
the average host star one) distribution in the Spitzer/IRAC 8 jim 
channel. The brightness distribution in Figure|3]is the mean of 
the brightness distribution trials accepted along the MCMC sim- 
ulation. Based on these brightness distribution trials, we also es- 
timate the dayside standard deviation to be about 0.025 in terms 
of relative brightness. 

In order to mitigate the model-dependence of our results, we 
use the three brightness models described in Subsection 2.2 All 
lead to the same conclusion: a brighter pattern to the North- 
East of the sub-stellar point and an eccentricity close to zero. 
In addition, our error bars are in agreement with our estimate of 
the brightness peak localization uncertainties reachable, about a 
dozen degrees in longitude and latitude, with current data qual- 
ity. We obtained this estimate in analysing simulated data. 

We hesitate to call the results a 2D map, because due to 
the limited data or data quality, we must use very simple mod- 
els which average the brightness distribution over large spatial 
scales. In other words, there is a degeneracy in the brightness 
anomaly distribution; the large feature in Figure[3] could in the 
future be spatially resolved into more complex structures. For 
the reason of limited data and geometry, the interpretation of 
Figure[3]has to focus on global trends, in particular: the presence 
of an asymmetrical hot spot. 

We present actually a time average large-scale brightness 
distribution of HD 189733b in Figure|3] left panel. In other 
words, Figure|3] shows the global pattern of HD 189733b bright- 
ness distribution based on eight snapshots taken from November 
2005 to June 2008. We focus this discussion on the HD 189733b 
dayside, as it actually results of the combination of both the 
phase curve and the occultations. HD 189733b dayside presents 
an offset hot spot. The eastward shifted distribution is in agree- 
ment with the literature: (1) with previous derivations, from 
HD 189733b phas e curve ([Knutson et al.||2007| ) and an eclipse 
timing constraint ( |Agol et al.||2010 ); and (2) atmospheric mod- 
els sugg esting a super-rotating equatorial jet (e.g.. Sho wman] 
et al.|2009j . In opposition, the north-south asymmetry is new. As 



we use large-scale, simple models to constrain the HD 189733b 
brightness distribution, the small-scale origin of this latitudi- 
nal asymmetry remains unconstrained. For example, a hot spot 
shifted to the East with a cold spot shifted to the South- West 
could mimic, within the current quality data, the photometry 
of our estimate of HD 189733b large-scale brightness distribu- 
tion. In addition, such a brightness distribution could a lso inter- 
estingly induce the specific feature observed in Knutson et al. 
( 2007 ) phase curve that is the phase curve extrema located be- 
tween the transit and the occultation. To explore further this pos- 
sibility and test the significance of more complex distribution 
we run additional MCMC simulations. These simulations use 
brightness models combining two of the simple brightness mod- 
els introduced in Subsec. [2!2| (i.e. two hot/cold zones). Their es- 
timate of HD 189733b large-scale brightness distribution are in 
excellent agreement with our results based on single "hot/cold 
spot" models. This confirms the low significance of small-scale 
pattern in current data and strengthen our previous results. 



Constraint on the HD 189733b orbital eccentricity 



Our work also yielded new system parameters in good agree- 
ment (to within less than l-cr) with previous work (e.g., Triaud 
et al.||2009l |Agol et al.]|2010| ). Notably, we refine the upper 
limit for the HD 189733b eccentricity: e < 0.0081 (95% con- 
fidence). Figure[3j right panel, presents the marginal posterior 
density distribution of e cos oj and e sin oj for our two analysis 
steps. The dashed lines refer to the first step using a uniform 
disk eclipse model, the solid lines to the second step based on 
the non-uniform brightness models. Figure[3j right panel, high- 
lights our degeneracy "lift" between the eccentricity (especially 
e cos <jj, see Sec. |2.2[ ) and the brightness distribution of an exo- 
planet. 



By a comparison with the uniform disk eclipse model, we 
see the "uniform timing offset", as a shift of e cos oj but no sig- 
nificant change in the e sin oj distribution. We note that the uni- 
form disk eclipse model estimation of e cos oj and e sin oj, while 
not used for the results , are in agreement with the literature (see 
e.g., [Triaud et al.|2009| ). 
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Fig. 3. Simultaneous estimation of HD 189733b global brightness distribution in IRAC 8 yum channel and eccentricity. Left: 
HD 189733b large-scale relative brightness distribution in IRAC 8 jim channel (estimated standard deviation in dayside: 0.025). 
Right: Marginal posterior densi ty pr obability distribution (68.3% and 95% confidence interva ls) o f ^coscl) and ^sinci; for the 
uniform brightness case (see Sec. |2.1[ dashed line) and the non-uniform brightness case (see Sec. 2^, solid line). 
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4. Summary and Discussion 

4.1. Methods and results 

We have performed a two-step analysis of the HD 189733b pub- 
lic eclipses obtained at 8 jim with Spitzer/IRAC. First, we have 
determined the system parameters using for occupation a uni- 
form disk eclipse model. This has led to the detection at the 
5cr level of occultation ingress/egress shapes that deviate from 
the ingress/egress of a uniform brightness disk. We have listed 
as major possible origins of these shapes, the exoplanet shape, 
brightness distribution or eccentricity. Based on the transits 
residuals we have rejected the exoplanet shape as possible ori- 
gin. Then, we have simultaneously analyzed our detection with 
HD 189733b phase curve from Knu tson et al.| ( |20Q7] ) to disentan- 
gle the contributions from the exoplanet brightness distribution 
in one hand and its eccentricity in the other hand. Our global 
analysis approach aims at mitigating the effects of systematics 
on individual analysis and at constraining consistently degenera- 
cies approached previously with strong assumptions. As a result, 
we have derived a highly precise upper limit on HD 189733b ec- 
centricity and have estimated its large-scale brightness distribu- 
tion in IRAC 8 fim channel. This one indicates a hot spot within 
the atmospheric layers probed as well as a double hemispheric 
asymmetry. As we have used large-scale, simple models to con- 
strain the HD 189733b brightness distribution, the small-scale 
origin of this latitudinal asymmetry remains unconstrained be- 
cause of the data quality and the planet orbit geometry. 



4.2. Comparison with \Majeau et a\.\\2012\ 



We obtain qualitatively similar results as Ml 2: an offset hot spot. 
Moreover, the estimations of the brightness distribution peak lo- 
cation are in good agreement (within less than Icr). We recall 
the main difference between both studies: in our work, 1) we 
propagate the system parameter uncertainties; 2) we do not con- 
strain the orbital eccentricity to zero but estimate it simultane- 
ously with the brightness distribution; 3) the starting point of 
our estimation of HD 189733b brightness distribution is its 5cr 
eclipse scanning. 

These approach differences lead to: 1) a significant differ- 
ence in the estimated uncertainty on the brightness peak lon- 
gitude, 21.8 ± 1.5°E compared to our 24 ± IFE; 2) differ- 
ent latitudinal offset significance, 3.1 ± 9.4°N compared to our 
17 ± 10°N. First, as introduced in Sec.[3]our error bars are con- 
sistent with our estimate of the brightness peak location uncer- 
tainty reachable with current data quality. Second, both M12's 
analysis and ours agree on the significant information of a lati- 
tudinal offset within the occultation ingress/egress: 13 ± 14.1°N 
(see Ml 2 Figures 1 and 3) compared to our 17 ± 10°N. However, 
in M12 the significance of this latitudinal offset decreases from 
13 ± 14.1°N to 3.1 ± 9.4°N when combined with the longitudinal 
information from the phase curve. 

4.3. Future 

Mapping distant worlds is a tremendous perspective that would 
dramatically improve our understanding of (exo)planets. As pre- 
viously introduced, observations of transiting exoplanets are 
sensitive to both their shape and their brightness distribution. In 
addition, the atmospheric layers actually scanned during occul- 
tation ingress/egress are wavelength-dependent. For that reason, 
with additional high quality occultation photometry at different 
wavelength we could improve our constraint on the brightness- 



eccentricity degeneracy and assess the stratified structure of an 
exoplanet atmosphere (e.g. higher wavelength, closer to the sub- 
stellar point and less extended hot spot). 

In addition the time variability of these spatial features could 
also be targeted. Indeed, as discussed in Rauscher et al.| ( [20Q7| ), 
/WST/NIRspec performance should allow high significance de- 
tection of modified occultation ingress/egress shapes based on a 
unique eclipse (for the grating centred at 4yL/m). The time vari- 
ability would then be asses from brightness derived from differ- 
ent occultations. 

In conclusion, the long-term challenge of constraining the 
exoplanet brightness distribution could ultimately lead to time- 
dependent three-dimensional map of distant worlds. 
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Table 2. AOR's description 



AORKEY(0') 



PI 



Publication Ref. 



Datasets^ (64x) 



16343552(0) 
20673792(P) 
22808832(0) 
22809088(0) 
22809344(0) 
22810112(0) 
24537600(0) 
27603456(0) 



D. Charbonneau 
D. Charbonneau 



E. Agol 



|Charbonneau e t al. ( 20081 
Knutson'eral., ( ,2007j 



Agol et aL| ( [20To] ) 



1359 
1319 



690 



Exposure time [s] Aperture [px] 



0.1 
0.4 



0.4 



3.6 
4.8 



4.8 



E. Agol 



Agol et al.| ( [20T0] ) 



690 



0.4 



22807296(T) 
22807552(T) 
22807808(T) 
24537856(T) 
276037 12(T) 
27773440(T) 

^ AORKEY target: T, O or P respectively transit, occultation or phase-curve. 

^ Present AOR are composed of datasets, each dataset corresponds to 64 individual subarray images of 32x32 px. 
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